function makeRunData_providerExits_ctsTime(runId)

addpath('../general_funcs');

config = getConfig_providerExits(runId);

dataId = config.dataId;

runDataFile           = config.runDataFile;
dynPredInputFile      = config.dynPredInputFile;

%%%%% Get covariates (including FEs)
covariates = config.covariates;

myX_tk_names   = covariates.myX_tk_names;
myX_tk_transfs = covariates.myX_tk_transfs;

myZ_names	= covariates.myZ_names;
myZ_transfs	= covariates.myZ_transfs;

myAnames	= covariates.myAnames;
myAtransfs	= covariates.myAtransfs;

myZ_geo_FE_Names  = covariates.myZ_geo_FE_Names;
time_FE_Names = covariates.time_FE_Names;

%%% Load files
data0 = loadDataFiles(dataId);
ctsTimeData  = data0.ctsTimeData;
periodicData = data0.periodicData;

%%%%%%%%%%%%%%%%%% LOAD STUFF FROM ctsTimeData %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NewConsumers        = ctsTimeData.NewConsumers;        % NumSpells x (K1*K2) - sparse
NewProviders         = ctsTimeData.NewProviders;         % NumSpells x K  - sparse
DeletedProviders     = ctsTimeData.DeletedProviders;     % NumSpells x K  - sparse
laggedConsumersOrig = ctsTimeData.laggedConsumersOrig; % NumSpells x K1 - sparse
laggedConsumersDest = ctsTimeData.laggedConsumersDest; % NumSpells x K2 - sparse
laggedProviders      = ctsTimeData.laggedProviders;      % NumSpells x K  - sparse

spellDurations    = ctsTimeData.spellDurations;    % NumSpells x 1
spell2period      = ctsTimeData.spell2period;      % NumSpells x 1 (integer between 1 and NumPeriods)

%%%%%%%%%%%%%%%%%% LOAD STUFF FROM periodicData %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
population        = periodicData.population;        % K x 1
surfaces          = periodicData.surfaces;          % K x 1
incomingCommuters = periodicData.incomingCommuters; % K x 1
touristBeds       = periodicData.touristBeds;       % K x 1
periodLabels = periodicData.periodLabels;
geoLabels   = periodicData.geoLabels;
A     = periodicData.A;
Z     = periodicData.Z;
X     = periodicData.X;
Anames     = periodicData.Anames;
Znames     = periodicData.Znames;
Xnames     = periodicData.Xnames;
surfaceName         = periodicData.surfaceName;
geoAggLevelNames    = periodicData.geoAggLevelNames;
geoAggLevelLabels   = periodicData.geoAggLevelLabels;
timeAggLevelNames   = periodicData.timeAggLevelNames;
timeAggLevelLabels = periodicData.timeAggLevelLabels;

% Read dimensions
NumSpells = size(NewProviders, 1);
NumPeriods = size(A,1);
K = size(Z,1);
K1 = K;
K2 = K;

%%% Make Y (=DeletedProviders)
Y = DeletedProviders;  % NumSpells x K (sparse)
clear DeletedProviders;


%%% Make logLambdaOffset (from pop at risk)
popAtRisk = laggedProviders;         % NumSpells x K2
logLambdaOffset = log(popAtRisk); % NumSpells x K2


% Reshape Y, installed bases and logLambdaOffset
Y                 = reshape(Y,                 [NumSpells*K  1]); % (NumSpells*K) x 1
laggedProviders      = reshape(laggedProviders,      [NumSpells*K2 1]); % (NumSpells*K2 x 1)
laggedConsumersOrig = reshape(laggedConsumersOrig, [NumSpells*K1 1]); % (NumSpells*K1 x 1)
laggedConsumersDest = reshape(laggedConsumersDest, [NumSpells*K2 1]); % (NumSpells*K2 x 1)
logLambdaOffset   = reshape(logLambdaOffset,   [NumSpells*K2, 1]); % (NumSpells*K) x 1

%%% Deal with cases where pop-at-risk is zero (then logLambdaOffset=-Inf, which causes issues)
% Avoid issues with logLambdaOffset=-Inf, while making sure that exp(logLambda) is equal to zero (computationallly speaking).
logLambdaOffset = max(logLambdaOffset, -2e4);

% Find corresponding tks that should be "deleted" and set the corresponding Y's to zero (it should be the case, but I need to make sure).
is_bad_tk  = laggedProviders <= 0; % (NumSpells*K) x 1
if any(Y(is_bad_tk) > 0); warning('Something strange here!'); end;
Y(is_bad_tk) = 0;



% Get myX_tk and myX_tk_names
[Xflag, X_idces] = ismember(myX_tk_names, Xnames);
if any(~Xflag); error('Could not find all Xnames.'); end;
if any(~strcmp(Xnames(X_idces), myX_tk_names)); error('Mistmatch'); end;
myX_tk         = X(:,:, X_idces);				% NumPeriods x K x NumX
myX_tk         = reshape(myX_tk, NumPeriods*K, size(myX_tk,3));	% (NumPeriods*K) x NumX
[myX_tk, myX_tk_names] = transformMatrix(myX_tk, myX_tk_transfs, myX_tk_names, reshape(repmat(surfaces', NumPeriods, 1), NumPeriods*K, 1),surfaceName);
	% (NumPeriods*K) x NumX and cell(1,NumX) of strings

% Get myZ and myZ_names
[Zflag, Z_idces] = ismember(myZ_names, Znames);
if any(~Zflag); error('Could not find all Znames.'); end;
if any(~strcmp(Znames(Z_idces), myZ_names)); error('Mistmatch'); end;
myZ         = Z(:, Z_idces);				% K x NumZ
[myZ, myZ_names] = transformMatrix(myZ, myZ_transfs, myZ_names, surfaces, surfaceName);
	% K x NumZ and cell(1,NumZ) of strings

% Get myA and myAnames
[Aflag, A_idces] = ismember(myAnames, Anames);
if any(~Aflag); error('Could not find all Anames.'); end;
if any(~strcmp(Anames(A_idces), myAnames)); error('Mistmatch'); end;
myA         = A(:, A_idces);				% NumPeriods x NumA
	% NumPeriods x NumA and cell(1,NumA) of strings

% Add intercept to myA
myA = [ones(NumPeriods,1) myA];
myAnames = [{'Intercept'} myAnames];
myAtransfs = [-1 myAtransfs];


%% Make time_FEs (I'm not using timeAggLevelNames and timeAggLevelLabels here but I could...)
[flag, idxes] = ismember(time_FE_Names, Anames);
if any(~flag); error('Could not find all time FE names.'); end;
time_FEs = A(:,idxes);
Num_time_FE_sets = length(time_FE_Names);
Num_time_FE_vals = zeros(1,Num_time_FE_sets);
time_FE_labels = cell(1,Num_time_FE_sets);
for aa = 1:Num_time_FE_sets
	time_FEs_aa = time_FEs(:,aa);
	time_FE_Name_aa = time_FE_Names{aa};
	
	unq_time_FEs_aa = sort(unique(time_FEs_aa));
	time_FE_labels_aa = strcat(time_FE_Name_aa, sprintfc('_is_%d', unq_time_FEs_aa));
	[flag,newTime_FEs_aa] = ismember(time_FEs_aa, unq_time_FEs_aa);
	
	time_FEs(:,aa) = newTime_FEs_aa;
	Num_time_FE_vals(aa) = length(unq_time_FEs_aa);
	time_FE_labels{aa} = time_FE_labels_aa;
end

%% Make  myZ_geo_FEs + get myZ_geo_FE_labels +  myZ_Num_geo_FE_vals
[flag, idxes] = ismember(myZ_geo_FE_Names, Znames);
if any(~flag); error('Could not find all geo FE names.'); end;
myZ_geo_FEs = Z(:,idxes);
Num_geo_FE_sets = length(myZ_geo_FE_Names);
myZ_Num_geo_FE_vals = zeros(1,Num_geo_FE_sets);
myZ_geo_FE_labels = cell(1,Num_geo_FE_sets);
for aa = 1:Num_geo_FE_sets
	geo_FE_Name_aa = myZ_geo_FE_Names{aa};
	[flag,geoidx] = ismember(geo_FE_Name_aa, geoAggLevelNames);
	if ~flag; error('Could not find geo_FE_name among geoAggLevelNames'); end;
	geo_FE_labels_aa = geoAggLevelLabels{geoidx};
	myZ_Num_geo_FE_vals(aa) = length(geo_FE_labels_aa);
	myZ_geo_FE_labels{aa} = geo_FE_labels_aa;
end

%%% Add treatment for initial installed base %%%
if config.initialBaseTreatment
	laggedProviders = reshape(laggedProviders, [NumSpells K]);  % NumSpells x K
	initialSeed = laggedProviders(1,:)'; % K x 1
	laggedProviders = reshape(laggedProviders, [NumSpells*K 1]); % (NumSpells*K) x 1
	[myZ, myZ_names] = addInitialBaseCovariate(initialSeed, myZ, myZ_names, 'Providers');
	clear initialSeed;
end

%%%% Make Xdynamic_staticInputs and Xparts_dynamic (Xparts that depend on installed base)
myX_tk = reshape(myX_tk, [NumPeriods, K, size(myX_tk,2)]);	% NumPeriods x K x NumX
myX_tk = myX_tk(spell2period,:,:); % NumSpells x K x NumX
myX_tk = reshape(myX_tk, [NumSpells*K, size(myX_tk,3)]);	% (NumSpells*K) x NumX

Xparts_static{1}.X      = myX_tk;             % (NumSpells*K) x NumXstatic
Xparts_static{1}.Xnames = myX_tk_names;       % cell(1,NumXstatic) of strings
Xparts_static{1}.X_FEs        = zeros(NumSpells*K,0); % (NumSpells*K) x NumX_FEs
Xparts_static{1}.X_FE_labels  = {};           % cell(1, NumX_FEs)
Xparts_static{1}.NumX_FE_vals = zeros(1,0);   % 1 x NumX_FEs

Xdynamic_staticInputs.Xparts_static = Xparts_static;
if config.laggedBaseTransf.code == 1 || config.laggedBaseTransf.code == 3
	Xdynamic_staticInputs.extraArgs     = {};
end
if config.laggedBaseTransf.code == 2 || config.laggedBaseTransf.code == 4  || config.laggedBaseTransf.code == 5
	if strcmp(config.laggedBaseTransf.denom, 'Area')
		surface_rep = reshape(repmat(surfaces', [NumSpells, 1]), [NumSpells*K 1]); % (NumSpells*K) x 1
		Xdynamic_staticInputs.extraArgs = {log(surface_rep)};
	end
end

installedBases.laggedProviders      = laggedProviders;
installedBases.laggedConsumersOrig = laggedConsumersOrig;
installedBases.laggedConsumersDest = laggedConsumersDest;
Xparts_dynamic = make_Xdynamic_providerExits(config, Xdynamic_staticInputs, installedBases, true);

%%%%% Make Xparts and dims %%%%%
dimType = 2;
Xparts = {};
dims1 = {};
Xpart_2_NumX = [];
Xpart_2_NumX_FEs = [];
Xpart_2_Num_FE_vals = {};
NumFEvals2Keep = {};

pidx = 0;

%%% T-K part:
pidx = pidx+1;
obj = Xparts_dynamic{1};
Xparts{pidx}              = obj;
dims1{pidx}               = NumSpells*K;
Xpart_2_NumX(pidx)        = size(obj.X,2);
Xpart_2_NumX_FEs(pidx)    = size(obj.X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = obj.NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(obj.NumX_FE_vals) - length(obj.NumX_FE_vals);

clear X obj;

%%% T part:
X      = myA(spell2period,:);            % NumSpells x NumX
Xnames = myAnames;                       % 1 x NumX
X_FEs        = time_FEs(spell2period,:); % NumSpells x NumX_FEs
X_FE_labels  = time_FE_labels;           % cell(1, NumX_FEs)
NumX_FE_vals = Num_time_FE_vals;         % 1 x NumX_FEs

pidx = pidx+1;
obj.X                     = X;
obj.Xnames                = Xnames;
obj.X_FEs                 = X_FEs;
obj.X_FE_labels           = X_FE_labels;
obj.NumX_FE_vals          = NumX_FE_vals;
Xparts{pidx}              = obj;
dims1{pidx}               = NumSpells;
Xpart_2_NumX(pidx)        = size(X,2);
Xpart_2_NumX_FEs(pidx)    = size(X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(NumX_FE_vals) - length(NumX_FE_vals);

clear X X_FEs NumX_FE_vals obj;

%%% K part:
X      = myZ;        % K x NumX
Xnames = myZ_names;  % 1 x NumX
X_FEs        = myZ_geo_FEs;  % K x NumX_FEs
X_FE_labels  = myZ_geo_FE_labels;   % cell(1, NumX_FEs)
NumX_FE_vals = myZ_Num_geo_FE_vals; % 1 x NumX_FEs

pidx = pidx+1;
obj.X                     = X;
obj.Xnames                = Xnames;
obj.X_FEs                 = X_FEs;
obj.X_FE_labels           = X_FE_labels;
obj.NumX_FE_vals          = NumX_FE_vals;
Xparts{pidx}              = obj;
dims1{pidx}               = K;
Xpart_2_NumX(pidx)        = size(X,2);
Xpart_2_NumX_FEs(pidx)    = size(X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(NumX_FE_vals) - length(NumX_FE_vals);

clear X X_FEs NumX_FE_vals obj;


%%%%% Make data object
% Calculate NumParams
NumParts = pidx;
NumParams = 0;
for ii = 1:NumParts
	NumParams = NumParams + size(Xparts{ii}.X,2);
	for ff = 1:size(Xparts{ii}.X_FEs,2)
		NumX_FE_vals_if = Xparts{ii}.NumX_FE_vals(ff);
		NumParams = NumParams + NumX_FE_vals_if - 1;
	end
end

% Make dims
dims.dimType  = dimType;
dims.NumParts = NumParts;
dims.NumObs = size(Y,1);
dims.NumParams      = NumParams;
dims.dims1               = dims1;
dims.Xpart_2_NumX        = Xpart_2_NumX;
dims.Xpart_2_NumX_FEs    = Xpart_2_NumX_FEs;
dims.Xpart_2_Num_FE_vals = Xpart_2_Num_FE_vals;
dims.NumFEvals2Keep      = NumFEvals2Keep;

% Make LL_offset
LL_offset = 0;

% Make data
clear data; % (Clear previous data object)
data.Y               = Y;                   % (NumSpells*K) x 1
data.logLambdaOffset = logLambdaOffset;     % (NumSpells*K) x 1
data.LL_offset       = LL_offset;           % scalar
data.Xparts          = Xparts;              % cell(3,1)
data.dims            = dims;                % cell(3,1)
data.periodLabels    = periodLabels;        % cell(NumPeriods,1)
data.geoLabels       = geoLabels;           % cell(K,1)
data.surfaces        = surfaces;            % K x 1
data.population      = population;          % K x 1
data.incomingCommuters = incomingCommuters; % K x 1
data.touristBeds       = touristBeds;       % K x 1
data.spellDurations  = spellDurations;      % NumSpells x 1
data.spell2period    = spell2period;        % NumSpells x 1 (values between 1 and NumPeriods)

% Save to output file
save(runDataFile, 'data', '-v7.3');

%%% Save things that are useful for dynamic prediction and counterfactuals
% Make Xdynamic_staticInputs at day level (not spell level)
period2firstSpell = accumarray(spell2period, [1:NumSpells]', [NumPeriods 1], @min); % NumPeriods x 1
if config.laggedBaseTransf.code == 2 || config.laggedBaseTransf.code == 3 || config.laggedBaseTransf.code == 4  || config.laggedBaseTransf.code == 5
	% Xparts_static{1}
	myX     = Xdynamic_staticInputs.Xparts_static{1}.X;     % (NumSpells*K) x NumXstatic
	myX_FEs = Xdynamic_staticInputs.Xparts_static{1}.X_FEs; % (NumSpells*K) x NumX_FEs
	
	myX     = reshape(myX,     [NumSpells K size(myX, 2)]);     % NumSpells x K x NumXstatic
	myX_FEs = reshape(myX_FEs, [NumSpells K size(myX_FEs, 2)]); % NumSpells x K x NumX_FEs
	
	myX     = myX(period2firstSpell,:,:);     % NumPeriods x K x NumXstatic
	myX_FEs = myX_FEs(period2firstSpell,:,:); % NumPeriods x K x NumX_FEs
	
	myX     = reshape(myX, [NumPeriods*K size(myX, 3)]);         % (NumPeriods*K) x NumXstatic
	myX_FEs = reshape(myX_FEs, [NumPeriods*K size(myX_FEs, 3)]); % (NumPeriods*K) x NumX_FEs
	
	Xdynamic_staticInputs.Xparts_static{1}.X     = myX;
	Xdynamic_staticInputs.Xparts_static{1}.X_FEs = myX_FEs;
	
	% extraArgs
	if config.laggedBaseTransf.code ~= 3
		tmp = Xdynamic_staticInputs.extraArgs{1}; % (NumSpells*K) x 1
		tmp = reshape(tmp, [NumSpells K]); % NumSpells x K
		tmp = tmp(period2firstSpell,:); % NumPeriods x K
		tmp = tmp(:); % (NumPeriods*K) x 1
		Xdynamic_staticInputs.extraArgs{1} = tmp; % (NumPeriods*K) x 1
	end
end

base.Y             = data.Y; % (sparse) NumSpells x K
base.Xdynamic_staticInputs = Xdynamic_staticInputs; % cell(NumParts,1)
base.laggedConsumersOrig = reshape(laggedConsumersOrig, [NumSpells K]);  % NumSpells x K
base.laggedConsumersDest = reshape(laggedConsumersDest, [NumSpells K]);  % NumSpells x K
base.laggedProviders      = reshape(laggedProviders,  [NumSpells K]); % NumSpells x K
base.popAtRisk     = popAtRisk; % NumSpells x K
save(dynPredInputFile, 'base', '-v7.3');

clear;
end